home *** CD-ROM | disk | FTP | other *** search
/ FishMarket 1.0 / FishMarket v1.0.iso / fishies / 376-400 / disk_386 / xlispstat / src2.lzh / XLisp-Stat / linalg.c < prev    next >
C/C++ Source or Header  |  1990-10-04  |  24KB  |  939 lines
  1. /* linalg - Lisp interface for basic linear algebra routines           */
  2. /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney                  */
  3. /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz    */
  4. /* You may give out copies of this software; for conditions see the    */
  5. /* file COPYING included with this distribution.                       */
  6.  
  7. #include <stdlib.h>
  8. #include "xlisp.h"
  9. #include "osdef.h"
  10. #ifdef ANSI
  11. #include "xlproto.h"
  12. #include "xlsproto.h"
  13. #else
  14. #include "xlfun.h"
  15. #include "xlsfun.h"
  16. #endif ANSI
  17. #include "xlsvar.h"
  18.  
  19. #ifdef ANSI
  20. LVAL vector_to_data(Vector,int,int,int),matrix_to_data(Matrix,int,int,int),
  21.      kernel(int),
  22.      add_contour_point(int,int,int,int,RVector,RVector,RMatrix,double,LVAL);
  23. char *allocate(unsigned,unsigned);
  24. void free_alloc(void *),baddata(LVAL),badmatrix(LVAL),badsquarematrix(LVAL),
  25.      badsequence(LVAL),
  26.      get_smoother_data(Vector *,Vector *,int *,Vector *,Vector *,int *,int);
  27. int data_mode(LVAL);
  28. Vector data_to_vector(LVAL,int);
  29. Matrix data_to_matrix(LVAL,int);
  30. #else
  31. LVAL vector_to_data(),matrix_to_data(),
  32.      kernel(),
  33.      add_contour_point();
  34. char *allocate();
  35. void free_alloc(),baddata(),badmatrix(),badsquarematrix(),
  36.      badsequence(),
  37.      get_smoother_data();
  38. int data_mode();
  39. Vector data_to_vector();
  40. Matrix data_to_matrix();
  41. #endif ANSI
  42.  
  43. /************************************************************************/
  44. /**                                                                    **/
  45. /**                 Storage Allocation Functions                       **/
  46. /**                                                                    **/
  47. /************************************************************************/
  48.  
  49. static char *allocate(n, m)
  50.     unsigned n, m;
  51. {
  52.   char *p = calloc(n, m);
  53.   if (p == nil) xlfail("allocation failed");
  54.   return(p);
  55. }
  56.  
  57. static void free_alloc(p)
  58. /*    char*/ void *p;/* changed JKL */
  59. {
  60.   if (p != nil) free(p);
  61. }
  62.  
  63. IVector ivector(n)
  64.     unsigned n;
  65. {
  66.   return((IVector) allocate(n, sizeof(int)));
  67. }
  68.  
  69. RVector rvector(n)
  70.     unsigned n;
  71. {
  72.   return((RVector) allocate(n, sizeof(double)));
  73. }
  74.  
  75. CVector cvector(n)
  76.     unsigned n;
  77. {
  78.   return((CVector) allocate(n, sizeof(Complex)));
  79. }
  80.  
  81. void free_vector(v) Vector v; { free_alloc(v); }
  82.  
  83. IMatrix imatrix(n, m)
  84.     unsigned n, m;
  85. {
  86.   int i;
  87.   IMatrix mat = (IMatrix) allocate(n, sizeof(IVector));
  88.   for (i = 0; i < n; i++) mat[i] = (IVector) allocate(m, sizeof(int));
  89.   return(mat);
  90. }
  91.  
  92. RMatrix rmatrix(n, m)
  93.     unsigned n, m;
  94. {
  95.   int i;
  96.   RMatrix mat = (RMatrix) allocate(n, sizeof(RVector));
  97.   for (i = 0; i < n; i++) mat[i] = (RVector) allocate(m, sizeof(double));
  98.   return(mat);
  99. }
  100.  
  101. CMatrix cmatrix(n, m)
  102.     unsigned n, m;
  103. {
  104.   int i;
  105.   CMatrix mat = (CMatrix) allocate(n, sizeof(CVector));
  106.   for (i = 0; i < n; i++) mat[i] = (CVector) allocate(m, sizeof(Complex));
  107.   return(mat);
  108. }
  109.  
  110. void free_matrix(mat, n)
  111.     Matrix mat;
  112.     int n;
  113. {
  114.   int i;
  115.   
  116.   if (mat != nil) for (i = 0; i < n; i++) free_alloc(mat[i]);
  117.   free_alloc(mat);
  118. }
  119.  
  120. /************************************************************************/
  121. /**                                                                    **/
  122. /**             Lisp to/from C Data Translation Functions              **/
  123. /**                                                                    **/
  124. /************************************************************************/
  125.  
  126. static void baddata(data)
  127.     LVAL data;
  128. {
  129.   xlerror("bad linear algebra data", data);
  130. }
  131.  
  132. static void badmatrix(data)
  133.     LVAL data;
  134. {
  135.   xlerror("not a matrix", data);
  136. }
  137.  
  138. static void badsquarematrix(data)
  139.     LVAL data;
  140. {
  141.   xlerror("not a square matrix", data);
  142. }
  143.  
  144. static void badsequence(data)
  145.     LVAL data;
  146. {
  147.   xlerror("not a sequence", data);
  148. }
  149.  
  150. static int data_mode(data)
  151.     LVAL data;
  152. {
  153.   LVAL item;
  154.   int i, n;
  155.   int mode = IN;
  156.   
  157.   if (! consp(data) && ! vectorp(data) && ! displacedarrayp(data))
  158.     baddata(data);
  159.     
  160.   if (! consp(data)) data = arraydata(data);
  161.   n = seqlen(data);
  162.   for (i = 0; i < n; i++) {
  163.     item = getnextelement(&data, i);
  164.     if (! realp(item) && ! complexp(item)) baddata(item);
  165.     if (floatp(item) && mode == IN) mode = RE;
  166.     else if (complexp(item)) mode = CX;
  167.   }
  168.   return(mode);
  169. }
  170.  
  171. static Vector data_to_vector(data, mode)
  172.     LVAL data;
  173.     int mode;
  174. {
  175.   LVAL item;
  176.   int i, n;
  177.   Vector v;
  178.   IVector iv;
  179.   RVector rv;
  180.   CVector cv;
  181.   
  182.   if (! sequencep(data)) baddata(data);
  183.   n = seqlen(data);
  184.   
  185.   switch (mode) {
  186.   case IN: iv = ivector(n); v = (Vector) iv; break;
  187.   case RE: rv = rvector(n); v = (Vector) rv; break;
  188.   case CX: cv = cvector(n); v = (Vector) cv; break;
  189.   }
  190.   
  191.   for (i = 0; i < n; i++) {
  192.     item = getnextelement(&data, i);
  193.     switch (mode) {
  194.     case IN: iv[i] = getfixnum(item);   break;
  195.     case RE: rv[i] = makedouble(item);  break;
  196.     case CX: cv[i] = makecomplex(item); break;
  197.     }
  198.   }
  199.   return(v);
  200. }
  201.  
  202. static Matrix data_to_matrix(data, mode)
  203.     LVAL data;
  204.     int mode;
  205. {
  206.   LVAL item;
  207.   int i, j, n, m;
  208.   Matrix mat;
  209.   IMatrix imat;
  210.   RMatrix rmat;
  211.   CMatrix cmat;
  212.   
  213.   if (! matrixp(data)) baddata(data);
  214.   n = numrows(data); m = numcols(data);
  215.   
  216.   switch (mode) {
  217.   case IN: imat = imatrix(n, m); mat = (Matrix) imat; break;
  218.   case RE: rmat = rmatrix(n, m); mat = (Matrix) rmat; break;
  219.   case CX: cmat = cmatrix(n, m); mat = (Matrix) cmat; break;
  220.   }
  221.   
  222.   data = arraydata(data);
  223.   for (i = 0; i < n; i++) {
  224.     for (j = 0; j < m; j++) {
  225.       item = getelement(data, i * m + j);
  226.       switch (mode) {
  227.       case IN: imat[i][j] = getfixnum(item);   break;
  228.       case RE: rmat[i][j] = makedouble(item);  break;
  229.       case CX: cmat[i][j] = makecomplex(item); break;
  230.       }
  231.     }
  232.   }
  233.   return(mat);
  234. }
  235.  
  236. static LVAL vector_to_data(v, n, mode, as_list)
  237.     Vector v;
  238.     int n, mode, as_list;
  239. {
  240.   LVAL data;
  241.   int i;
  242.   IVector iv = (IVector) v;
  243.   RVector rv = (RVector) v;
  244.   CVector cv = (CVector) v;
  245.   
  246.   xlsave1(data);
  247.   data = newvector(n);
  248.   for (i = 0; i < n; i++) {
  249.     switch (mode) {
  250.     case IN: setelement(data, i, cvfixnum((FIXTYPE) iv[i])); break;
  251.     case RE: setelement(data, i, cvflonum((FLOTYPE) rv[i])); break;
  252.     case CX: setelement(data, i, cvcomplex(cv[i]));          break;
  253.     }
  254.   }
  255.   if (as_list) data = coerce_to_list(data);
  256.   xlpop();
  257.   return(data);
  258. }
  259.  
  260. static LVAL matrix_to_data(mat, n, m, mode)
  261.     Matrix mat;
  262.     int n, m, mode;
  263. {
  264.   LVAL data, dim, data_matrix;
  265.   int i, j;
  266.   IMatrix imat = (IMatrix) mat;
  267.   RMatrix rmat = (RMatrix) mat;
  268.   CMatrix cmat = (CMatrix) mat;
  269.   
  270.   xlstkcheck(2);
  271.   xlsave(dim);
  272.   xlsave(data_matrix);
  273.   
  274.   dim = integer_list_2(n, m);
  275.   data_matrix = newarray(dim, NIL, NIL);
  276.   data = arraydata(data_matrix);
  277.   for (i = 0; i < n; i++) {
  278.     for (j = 0; j < m; j++) {
  279.       switch (mode) {
  280.       case IN: setelement(data, i * m + j, cvfixnum((FIXTYPE) imat[i][j])); break;
  281.       case RE: setelement(data, i * m + j, cvflonum((FLOTYPE) rmat[i][j])); break;
  282.       case CX: setelement(data, i * m + j, cvcomplex(cmat[i][j]));          break;
  283.       }
  284.     }
  285.   }
  286.   xlpopn(2);
  287.   return(data_matrix);
  288. }
  289.  
  290. /************************************************************************/
  291. /**                                                                    **/
  292. /**                  Machine Epsilon Determination                     **/
  293. /**                                                                    **/
  294. /************************************************************************/
  295.  
  296. double macheps()
  297. {
  298.   static int calculated = FALSE;
  299.   static double epsilon = 1.0;
  300.   
  301.   if (! calculated)
  302.     while (1.0 + epsilon / 2.0 != 1.0) epsilon = epsilon / 2.0;
  303.   calculated = TRUE;
  304.   return(epsilon);
  305. }
  306.  
  307. /************************************************************************/
  308. /**                                                                    **/
  309. /**            Lisp Interfaces to Linear Algebra Routines              **/
  310. /**                                                                    **/
  311. /************************************************************************/
  312.  
  313. LVAL xslu_decomp()
  314. {
  315.   LVAL data, result, temp;
  316.   Matrix mat;
  317.   IVector iv;
  318.   double d;
  319.   int n, m, mode, singular;
  320.   
  321.   data = xlgetarg();
  322.   xllastarg();
  323.   
  324.   if (! matrixp(data)) badmatrix(data);
  325.   n = numrows(data);
  326.   m = numcols(data);
  327.   if (n != m || n <= 0 || m <= 0) badsquarematrix(data);
  328.   
  329.   xlsave1(result);
  330.   result = mklist(4, NIL); /* redundant assignment to nil in macro JKL */
  331.   temp = result;
  332.   
  333.   mode = data_mode(data);
  334.   if (mode == IN) mode = RE;
  335.   mat = data_to_matrix(data, mode);
  336.   iv = ivector(n);
  337.   singular = crludcmp(mat, n, iv, mode, &d);
  338.   rplaca(temp, matrix_to_data(mat, n, n, mode));    temp = cdr(temp);
  339.   rplaca(temp, vector_to_data((Vector)iv, n, IN, FALSE));   temp = cdr(temp);
  340.   rplaca(temp, cvflonum((FLOTYPE) d));              temp = cdr(temp);
  341.   rplaca(temp, (singular) ? s_true : NIL);
  342.  
  343.   free_matrix(mat, n);
  344.   free_vector((Vector)iv); /* casts added JKL */
  345.   xlpop();
  346.   return(result);
  347. }
  348.  
  349. LVAL xslu_solve()
  350. {
  351.   LVAL ludecomp, la, lindx, lb, result;
  352.   Matrix a;
  353.   Vector b, indx;
  354.   int n, m, a_mode, b_mode, mode, singular;
  355.   
  356.   ludecomp = xlgalist();
  357.   lb = xlgetarg();
  358.   xllastarg();
  359.   
  360.   la = (consp(ludecomp)) ? car(ludecomp) : NIL;
  361.   lindx = (consp(cdr(ludecomp))) ? car(cdr(ludecomp)) : NIL;
  362.   if (! matrixp(la)) badmatrix(la);
  363.   n = numrows(la);
  364.   m = numcols(la);
  365.   
  366.   if (n != m || n <= 0 || m <= 0) badsquarematrix(la);
  367.   
  368.   if (! sequencep(lindx)) badsequence(lindx);
  369.   if (data_mode(lindx) != IN) xlerror("not an integer sequence", lindx);
  370.   if (! sequencep(lb)) badsequence(lb);
  371.   if (seqlen(lb) != n) xlerror("bad sequence length", lb);
  372.   
  373.   a_mode = data_mode(la);
  374.   b_mode = data_mode(lb);
  375.   mode = max(a_mode, b_mode);
  376.   if (mode == IN) mode = RE;
  377.   
  378.   a = data_to_matrix(la, mode);
  379.   indx = data_to_vector(lindx, IN);
  380.   b = data_to_vector(lb, mode);
  381.   
  382.   singular = crlubksb(a, n, (IVector)indx, b, mode);/* cast added JKL */
  383.   result = vector_to_data(b, n, mode, listp(lb));
  384.   free_matrix(a, n);
  385.   free_vector(indx);
  386.   free_vector(b);
  387.   
  388.   if (singular) xlfail("matrix is (numerically) singular");
  389.   return(result);
  390. }
  391.  
  392. LVAL xslu_determinant()
  393. {
  394.   LVAL data, result;
  395.   Matrix mat;
  396.   CMatrix cmat;
  397.   RMatrix rmat;
  398.   IVector iv;
  399.   double d, rd1, d2, magn;
  400.   Complex cd1;
  401.   int m, n, i, mode;
  402.   
  403.   data = xlgetarg();
  404.   xllastarg();
  405.   
  406.   if (! matrixp(data)) badmatrix(data);
  407.   m = numrows(data);
  408.   n = numcols(data);
  409.   if (n != m || n <= 0 || m <= 0) badsquarematrix(data);
  410.     
  411.   mode = data_mode(data);
  412.   if (mode == IN) mode = RE;
  413.   mat = data_to_matrix(data, mode);
  414.   iv = ivector(n);
  415.   crludcmp(mat, n, iv, mode, &d);
  416.  
  417.   switch (mode) {
  418.   case RE: 
  419.     rmat = (RMatrix) mat;
  420.     rd1 = d;
  421.     d2 = 0.0;
  422.     for (i = 0; i < n; i++) {
  423.       if ((magn = fabs(rmat[i][i])) == 0.0) {
  424.         rd1 = 0.0;
  425.         break;
  426.       }
  427.       rd1 *= rmat[i][i] / magn;
  428.       d2 += log(magn);
  429.     }
  430.     result = cvflonum((FLOTYPE) rd1 * exp(d2));
  431.     break;
  432.   case CX:
  433.     cmat = (CMatrix) mat;
  434.     cd1 = cart2complex(d, 0.0);
  435.     d2 = 0.0;
  436.     for (i = 0; i < n; i++) {
  437.       if ((magn = modulus(cmat[i][i])) == 0.0) {
  438.         cd1 = cart2complex(0.0, 0.0);
  439.         break;
  440.       }
  441.       cd1 = polar2complex(modulus(cd1), phase(cd1) + phase(cmat[i][i]));
  442.       d2 += log(magn);
  443.     }
  444.     result = cvcomplex(cmul(cd1, cart2complex(exp(d2), 0.0)));
  445.     break;
  446.   }
  447.   
  448.   free_matrix(mat, n);
  449.   free_vector((Vector)iv);/* cast added JKL */
  450.   return(result);
  451. }
  452.  
  453. LVAL xslu_inverse()
  454. {
  455.   LVAL data, result;
  456.   Matrix mat, inv;
  457.   CMatrix cinv;
  458.   RMatrix rinv;
  459.   CVector cv;
  460.   RVector rv;
  461.   Vector v;
  462.   IVector iv;
  463.   double d;
  464.   int m, n, i, j, mode, singular;
  465.   
  466.   data = xlgetarg();
  467.   xllastarg();
  468.   
  469.   if (! matrixp(data)) badmatrix(data);
  470.   m = numrows(data);
  471.   n = numcols(data);
  472.   if (n != m || n <= 0 || m <= 0) badsquarematrix(data);
  473.     
  474.   mode = data_mode(data);
  475.   if (mode == IN) mode = RE;
  476.   mat = data_to_matrix(data, mode);
  477.   iv = ivector(n);
  478.   inv = (mode == RE) ? (Matrix) rmatrix(n,n) : (Matrix) cmatrix(n,n);
  479.   rinv = (RMatrix) inv;
  480.   cinv = (CMatrix) inv;
  481.   v = (mode == RE) ? (Vector) rvector(n) : (Vector) cvector(n);
  482.   rv = (RVector) v;
  483.   cv = (CVector) v;
  484.   
  485.   singular = crludcmp(mat, n, iv, mode, &d);
  486.  
  487.   if (! singular) {
  488.     for (j = 0; j < n; j++) {
  489.       for (i = 0; i < n; i++) {
  490.         if (mode == RE) rv[i] = 0.0;
  491.         else cv[i] = cart2complex(0.0, 0.0);
  492.       }
  493.       if (mode == RE) rv[j] = 1.0;
  494.       else cv[j] = cart2complex(1.0, 0.0);
  495.       
  496.       singular = singular || crlubksb(mat, n, iv, v, mode);
  497.       
  498.       for (i = 0; i < n; i++) {
  499.         if (mode == RE) rinv[i][j] = rv[i];
  500.         else cinv[i][j] = cv[i];
  501.       }
  502.     }
  503.     result = matrix_to_data(inv, n, n, mode);
  504.   }
  505.   
  506.   free_matrix(mat, n);
  507.   free_matrix(inv, n);
  508.   free_vector((Vector)iv);/* cast added JKL */
  509.   free_vector(v);
  510.   if (singular) xlfail("matrix is (numerically) singular");
  511.   return(result);
  512. }
  513.  
  514. LVAL xssv_decomp()
  515. {
  516.   LVAL data, result, temp;
  517.   Matrix mat, v;
  518.   Vector w;
  519.   int n, m, mode, converged;
  520.   
  521.   data = xlgetarg();
  522.   xllastarg();
  523.   
  524.   if (! matrixp(data)) badmatrix(data);
  525.   m = numrows(data);
  526.   n = numcols(data);
  527.   if (n <= 0 || m <= 0) baddata(data);
  528.   if (m < n) xlfail("number of rows less than number of columns");
  529.   
  530.   xlsave1(result);
  531.   result = mklist(4, NIL);
  532.   temp = result;
  533.   
  534.   mode = data_mode(data);
  535.   if (mode == IN) mode = RE;
  536.   if (mode == CX) xlfail("complex SVD not available yet");
  537.   
  538.   mat = data_to_matrix(data, mode);
  539.   w = (Vector) rvector(n);
  540.   v = (Matrix) rmatrix(n, n);
  541.                      /* casts added JKL */
  542.   converged = svdcmp((RMatrix)mat, m, n, (RVector)w, (RMatrix)v);
  543.   rplaca(temp, matrix_to_data(mat, m, n, mode));     temp = cdr(temp);
  544.   rplaca(temp, vector_to_data(w, n, mode, FALSE));   temp = cdr(temp);
  545.   rplaca(temp, matrix_to_data(v, n, n, mode));       temp = cdr(temp);
  546.   rplaca(temp, (converged) ? s_true : NIL);
  547.  
  548.   free_matrix(mat, m);
  549.   free_matrix(v, n);
  550.   free_vector(w);
  551.   xlpop();
  552.   return(result);
  553. }
  554.  
  555. LVAL xsqr_decomp()
  556. {
  557.   LVAL data, result, temp;
  558.   Matrix mat, v;
  559.   Vector jpvt;
  560.   int n, m, mode, pivot;
  561.   
  562.   data = xlgetarg();
  563.   pivot = (moreargs()) ? (xlgetarg() != NIL) : FALSE;
  564.   xllastarg();
  565.   
  566.   if (! matrixp(data)) badmatrix(data);
  567.   m = numrows(data);
  568.   n = numcols(data);
  569.   if (n <= 0 || m <= 0) baddata(data);
  570.   if (m < n) xlfail("number of rows less than number of columns");
  571.   
  572.   xlsave1(result);
  573.   result = mklist((pivot) ? 3 : 2, NIL);
  574.   temp = result;
  575.   
  576.   mode = data_mode(data);
  577.   if (mode == IN) mode = RE;
  578.   if (mode == CX) xlfail("complex QR decomposition not available yet");
  579.   
  580.   mat = data_to_matrix(data, mode);
  581.   v = (Matrix) rmatrix(n, n);
  582.   jpvt = (Vector) ivector(n);
  583.            /* casts added JKL */
  584.   qrdecomp((RMatrix)mat, m, n, (RMatrix)v, (IVector)jpvt, pivot);
  585.   rplaca(temp, matrix_to_data(mat, m, n, mode)); temp = cdr(temp);
  586.   rplaca(temp, matrix_to_data(v, n, n, mode)); temp = cdr(temp);
  587.   if (pivot) rplaca(temp, vector_to_data((Vector)jpvt, n, IN, TRUE));
  588.  
  589.   free_matrix((Matrix)mat, m);/* casts added JKL */
  590.   free_matrix((Matrix)v, n);
  591.   free_vector((Vector)jpvt);
  592.  
  593.   xlpop();
  594.   return(result);
  595. }
  596.  
  597. LVAL xschol_decomp()
  598. {
  599.   LVAL data, result;
  600.   Matrix mat;
  601.   int n, m, mode;
  602.   double maxadd, maxoffl;
  603.   
  604.   data = xlgetarg();
  605.   if (moreargs()) maxoffl = makedouble(xlgetarg());
  606.   else maxoffl = 0.0;
  607.   xllastarg();
  608.   
  609.   if (! matrixp(data)) badmatrix(data);
  610.   m = numrows(data);
  611.   n = numcols(data);
  612.   if (n != m || n <= 0 || m <= 0) badsquarematrix(data);
  613.     
  614.   mode = data_mode(data);
  615.   if (mode == IN) mode = RE;
  616.   if (mode == CX) xlfail("complex Cholesky not available yet");
  617.   
  618.   mat = data_to_matrix(data, mode);
  619.  
  620.   choldecomp((RMatrix)mat, n, maxoffl, &maxadd);/* cast added JKL */
  621.   result = mklist(2, NIL);
  622.   rplaca(result, matrix_to_data(mat, n, n, mode));
  623.   rplaca(cdr(result), cvflonum((FLOTYPE) maxadd));
  624.  
  625.   free_matrix(mat, m);
  626.   return(result);
  627. }
  628.  
  629. LVAL xsrcondest()
  630. {
  631.   LVAL data, result;
  632.   Matrix mat;
  633.   int n, m, mode;
  634.   double est;
  635.   
  636.   data = xlgetarg();
  637.   xllastarg();
  638.   
  639.   if (! matrixp(data)) badmatrix(data);
  640.   m = numrows(data);
  641.   n = numcols(data);
  642.   if (n != m || n <= 0 || m <= 0) badsquarematrix(data);
  643.     
  644.   mode = data_mode(data);
  645.   if (mode == IN) mode = RE;
  646.   if (mode == CX) xlfail("complex condition estimate not available yet");
  647.   
  648.   mat = data_to_matrix(data, mode);
  649.  
  650.   est = rcondest((RMatrix)mat, n);/* cast added JKL */
  651.   result = cvflonum((FLOTYPE) est);
  652.  
  653.   free_matrix(mat, m);
  654.   return(result);
  655. }
  656.  
  657. LVAL xsmake_rotation()
  658. {
  659.   LVAL s1, s2, result;
  660.   Vector x, y;
  661.   Matrix rot;
  662.   double alpha;
  663.   int n, use_alpha = FALSE;
  664.   
  665.   s1 = xsgetsequence();
  666.   s2 = xsgetsequence();
  667.   if (moreargs()) {
  668.     use_alpha = TRUE;
  669.     alpha = makedouble(xlgetarg());
  670.   }
  671.   xllastarg();
  672.   
  673.   n = seqlen(s1);
  674.   if (seqlen(s2) != n) xlfail("sequences not the same length");
  675.   
  676.   if (data_mode(s1) == CX || data_mode(s2) == CX)
  677.     xlfail("complex data not supported yet");
  678.   
  679.   x = data_to_vector(s1, RE);
  680.   y = data_to_vector(s2, RE);
  681.   
  682.   rot = (Matrix) rmatrix(n,n);/* casts added JKL */
  683.   make_rotation(n, (RMatrix)rot, (RVector)x, (RVector)y, use_alpha, alpha);
  684.   result = matrix_to_data(rot, n, n, RE);
  685.   
  686.   free_vector(x);
  687.   free_vector(y);
  688.   free_matrix(rot, n);
  689.   return(result);
  690. }
  691.  
  692. #define NS_DEFAULT 30
  693.  
  694. static void get_smoother_data(px, py, pn, pxs, pys, pns, is_reg)
  695.      Vector *px, *py, *pxs, *pys;
  696.      int *pn, *pns, is_reg;
  697. {
  698.   LVAL s1, s2, arg, sk_xvals = xlenter(":XVALS");
  699.   int n, ns, i, supplied; 
  700.   double xmin, xmax, *x, *xs;
  701.  
  702.   s1 = xsgetsequence();
  703.   if (is_reg) s2 = xsgetsequence();
  704.  
  705.   if (xlgetkeyarg(sk_xvals, &arg)) {
  706.     if (! sequencep(arg) && ! fixp(arg)) xlbadtype(arg);
  707.   }
  708.   else arg = NIL;
  709.  
  710.   ns = (fixp(arg)) ? getfixnum(arg) : seqlen(arg);
  711.   if (ns < 2) ns = NS_DEFAULT;
  712.   supplied = (sequencep(arg) && arg != NIL) ? TRUE : FALSE;
  713.  
  714.   n = seqlen(s1);
  715.   if (n <= 0) xlfail("sequence too short");
  716.   if (is_reg && seqlen(s2) != n) xlfail("sequences not the same length");
  717.   if (supplied && data_mode(arg) == CX) xlfail("data must be real");
  718.   if (data_mode(s1) == CX || (is_reg && data_mode(s2) == CX))
  719.     xlfail("data must be real");
  720.   
  721.   *px = data_to_vector(s1, RE);
  722.   *py = (is_reg) ? data_to_vector(s2, RE) : nil;
  723.   *pxs = (supplied) ? data_to_vector(arg, RE) : (Vector) rvector(ns);
  724.   *pys = (Vector) rvector(ns);
  725.   *pn = n;
  726.   *pns = ns;
  727.  
  728.   if (! supplied) {
  729.     x = (double *) *px;
  730.     xs = (double *) *pxs;
  731.     for (xmax = xmin = x[0], i = 1; i < *pn; i++) {
  732.       if (x[i] > xmax) xmax = x[i];
  733.       if (x[i] < xmin) xmin = x[i];
  734.     }
  735.     for (i = 0; i < *pns; i++)
  736.       xs[i] = xmin + (xmax - xmin) * ((double) i) / ((double) (*pns - 1));
  737.   }
  738. }
  739.  
  740. LVAL xsspline()
  741. {
  742.   LVAL result;
  743.   Vector x, y, work, xs, ys;
  744.   int n, ns, error;
  745.  
  746.   get_smoother_data(&x, &y, &n, &xs, &ys, &ns, TRUE);
  747.  
  748.   work = (Vector) rvector(2 * n);
  749.                         /* casts added JKL */
  750.   error = fit_spline(n, (RVector)x, (RVector)y, ns, (RVector)xs,
  751.                      (RVector)ys, (RVector)work);
  752.   xlsave1(result);
  753.   result = mklist(2, NIL);
  754.   rplaca(result, vector_to_data(xs, ns, RE, TRUE));
  755.   rplaca(cdr(result), vector_to_data(ys, ns, RE, TRUE));
  756.   xlpop();
  757.  
  758.   free_vector(x);
  759.   free_vector(y);
  760.   free_vector(xs);
  761.   free_vector(ys);
  762.   free_vector(work);
  763.  
  764.   if (error) xlfail("bad data for splines");
  765.   return(result);
  766. }
  767.  
  768. static LVAL kernel(is_reg)
  769.      int is_reg;
  770. {
  771.   LVAL warg, targ, result;
  772.   Vector x, y, xs, ys/*, wts, wds*/;
  773.   int n, ns, error, ktype;
  774.   double width;
  775.  
  776.   get_smoother_data(&x, &y, &n, &xs, &ys, &ns, is_reg);
  777.   if (! xlgetkeyarg(sk_width, &warg)) warg = NIL;
  778.   if (! xlgetkeyarg(sk_type, &targ)) targ = NIL;
  779.  
  780.   width = (fixp(warg) || floatp(warg)) ? makedouble(warg) : -1.0;
  781.   /*wts = (Vector) nil;
  782.     wds = (Vector) nil; not used JKL */
  783.   if (targ == xlenter("U")) ktype = 'U';
  784.   else if (targ == xlenter("T")) ktype = 'T';
  785.   else if (targ == xlenter("G")) ktype = 'G';
  786.   else ktype = 'B';
  787.                             /* casts added JKL */
  788.   error = kernel_smooth((RVector)x, (RVector)y, n, width, nil, nil,
  789.                         (RVector)xs, (RVector)ys, ns, ktype);
  790.   xlsave1(result); /* redundant equation of result to NIL in macro JKL */
  791.   result = mklist(2, NIL);
  792.   rplaca(result, vector_to_data(xs, ns, RE, TRUE));
  793.   rplaca(cdr(result), vector_to_data(ys, ns, RE, TRUE));
  794.   xlpop();
  795.  
  796.   free_vector(x);
  797.   free_vector(y);
  798.   free_vector(xs);
  799.   free_vector(ys);
  800.  
  801.   if (error) xlfail("bad data for splines");
  802.   return(result);
  803. }
  804.  
  805. LVAL xskernel_smooth() { return(kernel(TRUE));  }
  806. LVAL xskernel_dens()   { return(kernel(FALSE)); }
  807.  
  808. LVAL xsbase_lowess() 
  809. {
  810.   LVAL s1, s2, result;
  811.   int n, nsteps, error;
  812.   double f, delta;
  813.   Vector x, y, ys, rw, res;
  814.  
  815.   s1 = xsgetsequence();
  816.   s2 = xsgetsequence();
  817.   f = makedouble(xlgetarg());
  818.   nsteps = getfixnum(xlgafixnum());
  819.   delta = makedouble(xlgetarg());
  820.   xllastarg();
  821.  
  822.   n = seqlen(s1);
  823.   if (n <= 0) xlfail("sequence too short");
  824.   if (seqlen(s2) != n) xlfail("sequences not the same length");
  825.   if (data_mode(s1) == CX || data_mode(s2) == CX) xlfail("data must be real");
  826.   x = data_to_vector(s1, RE);
  827.   y = data_to_vector(s2, RE); 
  828.   ys = (Vector) rvector(n);
  829.   rw = (Vector) rvector(n);
  830.   res = (Vector) rvector(n);
  831.                  /* casts added JKL */
  832.   error = lowess((RVector)x, (RVector)y, n, f, nsteps, delta, (RVector)ys,
  833.                  (RVector)rw, (RVector)res);
  834.   result = vector_to_data(ys, n, RE, TRUE);
  835.  
  836.   free_vector(x);
  837.   free_vector(y);
  838.   free_vector(ys);
  839.   free_vector(rw);
  840.   free_vector(res);
  841.  
  842.   if (error) xlfail("bad data for splines");
  843.   return(result);
  844. }
  845.  
  846. static LVAL add_contour_point(i, j, k, l, x, y, z, v, result)
  847.     int i, j, k, l;
  848.     RVector x, y;
  849.     RMatrix z;
  850.     double v;
  851.     LVAL result;
  852. {
  853.   LVAL pt;
  854.   double p, q;
  855.   
  856.   if ((z[i][j] <= v && v < z[k][l]) || (z[k][l] <= v && v < z[i][j])) {
  857.     xlsave(pt);
  858.     pt = mklist(2, NIL);
  859.     p = (v - z[i][j]) / (z[k][l] - z[i][j]);
  860.     q = 1.0 - p;
  861.     rplaca(pt, cvflonum((FLOTYPE) (q * x[i] + p * x[k])));
  862.     rplaca(cdr(pt), cvflonum((FLOTYPE) (q * y[j] + p * y[l])));
  863.     result = cons(pt, result);
  864.     xlpop();
  865.   }
  866.   return(result);
  867. }
  868.  
  869. LVAL xssurface_contour()
  870. {
  871.   LVAL s1, s2, mat, result;
  872.   RVector x, y;
  873.   RMatrix z;
  874.   double v;
  875.   int i, j, n, m;
  876.   
  877.   s1 = xsgetsequence();
  878.   s2 = xsgetsequence();
  879.   mat = xsgetmatrix();
  880.   v = makedouble(xlgetarg());
  881.   xllastarg();
  882.     
  883.   n = seqlen(s1); m = seqlen(s2);
  884.   if (n != numrows(mat) || m != numcols(mat)) xlfail("dimensions do not match");
  885.   if (data_mode(s1) == CX || data_mode(s2) == CX || data_mode(mat) == CX)
  886.     xlfail("data must be real");
  887.   
  888.   x = (RVector) data_to_vector(s1, RE);
  889.   y = (RVector) data_to_vector(s2, RE);
  890.   z = (RMatrix) data_to_matrix(mat, RE);
  891.   
  892.   xlsave1(result);
  893.   result = NIL;
  894.   for (i = 0; i < n - 1; i++) {
  895.       for (j = 0; j < m - 1; j++) {
  896.         result = add_contour_point(i, j, i, j+1, x, y, z, v, result);
  897.         result = add_contour_point(i, j+1, i+1, j+1, x, y, z, v, result);
  898.         result = add_contour_point(i+1, j+1, i+1, j, x, y, z, v, result);
  899.         result = add_contour_point(i+1, j, i, j, x, y, z, v, result);
  900.       }
  901.   }
  902.   xlpop();
  903.   
  904.   free_vector((Vector)x);/* casts added JKL */
  905.   free_vector((Vector)y);
  906.   free_matrix((Matrix)z, n);
  907.   
  908.   return(result);
  909. }
  910.  
  911. LVAL xsfft()
  912. {
  913.   LVAL data, result;
  914.   CVector x;
  915.   RVector work;
  916.   int n, isign, as_list;
  917.   
  918.   data = xsgetsequence();
  919.   isign = (moreargs() && xlgetarg() != NIL) ? -1.0 : 1.0; 
  920.   xllastarg();
  921.   
  922.   /* check and convert the data */
  923.   n = seqlen(data);
  924.   if (n <= 0) xlfail("not enough data");
  925.   data_mode(data); /* checks that data are numbers */
  926.   as_list = (listp(data)) ? TRUE : FALSE;
  927.   x = (CVector) data_to_vector(data, CX);
  928.   work = rvector(4 * n + 15);
  929.  
  930.   cfft(n, (double *)x, (double *)work, isign); /* casts added JKL */
  931.  
  932.   /* free internal data and return result */
  933.   result = vector_to_data((Vector)x, n, CX, as_list);/* casts added JKL */
  934.   free_vector((Vector)x);
  935.   free_vector((Vector)work);
  936.   
  937.   return(result);
  938. }
  939.